Author(s): Olli Väinölä olli.vainola@student.oulu.fi
In this notebook is an small tutorial, how to create a Von Mises material with hardening. Equations are formulated into rate depended form. Code may or may not include some bugs..
Stress:
$\sigma = C : \epsilon^e$
$C$ is material tensor and $\epsilon$ total strain. Total strain is divided into elastic and plastic part:
$\epsilon = \epsilon^e + \epsilon^p$
now let's define a strain rate, which is strain increment divide with time increment $dt$
$\frac{d\epsilon}{dt} = \dot \epsilon = \dot \epsilon^e + \dot \epsilon^p$
Now same procedure for stress and substitute $\dot \epsilon^e$
$\dot \sigma = C : \dot \epsilon^e = C : (\dot \epsilon - \dot \epsilon^p )$
Only thing to do is to define yield function. Now we're using Von Mises material:
$f(\sigma - X, R(\alpha)) = \sqrt{3J_2(\sigma-X))} - R(\alpha)$ = 0
$f(\sigma, \sigma_y) = \sqrt{3J_2(\sigma))} - \sigma_y$ = 0
$J_2 = \frac{1}{2}s : s$
$s = \sigma - \frac{1}{3}\sigma I$
$I = eye(3)$
In [1]:
# imports
using PyPlot
using ForwardDiff
using NLsolve
Let's create a isotropic Hooke material.
In [7]:
"""
Create a isotropic Hooke material matrix C
More information: # http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm
Parameters
----------
E: Float
Elastic modulus
ν: Float
Poisson constant
Returns
-------
Array{Float64, (6,6)}
"""
function hookeStiffnessTensor(E, ν)
a = 1 - ν
b = 1 - 2*ν
c = 1 + ν
multiplier = E / (b * c)
return Float64[a ν ν 0 0 0;
ν a ν 0 0 0;
ν ν a 0 0 0;
0 0 0 b 0 0;
0 0 0 0 b 0;
0 0 0 0 0 b].*multiplier
end
# Pick material values
E = 200.0e3
ν = 0.3
C = hookeStiffnessTensor(E, ν)
Out[7]:
In [178]:
# using vectors with double contradiction
# http://www-2.unipv.it/compmech/teaching/available/const_mod/const_mod_mat-review_notation.pdf
M = [1 0 0 0 0 0;
0 1 0 0 0 0;
0 0 1 0 0 0;
0 0 0 2 0 0;
0 0 0 0 2 0;
0 0 0 0 0 2;]
"""
Equivalent tensile stress.
More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion
Section: Reduced von Mises equation for different stress conditions
Parameters
----------
σ: Array{Float64, 6}
Stress in Voigt notation
Returns
-------
Float
"""
function σₑ(σ)
s = σ[1:6] - 1/3 * sum([σ[1], σ[2], σ[3]]) * [1 1 1 0 0 0]'
return sqrt(3/2 * s' * M * s)[1]
end
# Some extra data for testing purposes ...
#ss = σ[1]
#return sqrt(sum(ss.^2))
"""
Von Mises Yield criterion
More info can be found from: http://csm.mech.utah.edu/content/wp-content/uploads/2011/10/9tutorialOnJ2Plasticity.pdf
Parameters
----------
σ: Array{Float64, 6}
Stress in Voigt notation
k: Float64
Material constant, Yield limit
Returns
-------
Float
"""
function vonMisesYield(σ, k)
σₑ(σ) - k
end
"""
Function for NLsolve. Inside this function are the functions where we want to find root.
Ψ is the yield function below. Functions defined here:
dσ - C (dϵ - dλ*dΨ/dσ) = 0
σₑ(σ) - k = 0
Parameters
----------
params: Array{Float64, 7}
Array containing values from solver
dϵ: Array{Float64, 6}
Strain rate vector in Voigt notation
C: Array{Float64, (6, 6)}
Material tensor
k: Float
Material constant, yield limit
Δt: Float
time increment
σ_begin:Array{Float64, 6}
Stress vector in Voigt notation
p: Float
Accumulated plastic strain
X: Array{Float64, 6}
Kinematic hardening tensor
R: Function
Calculates isotropic hardening as a function of accumulated plastic strain
α:Float
Accumulated kinematic evolution variable
dα: Function
Calculates kinematic evolution variables rate as a function of accumulated plastic slip rate
fX: Function
Calculates the kinematic
Returns
-------
Array{Float64, 7}, return values for solver
"""
function G(params, dϵ, C, σ_begin, p, X, R, fX, ϵp)
# Initializing wrapper
yield(pars) = vonMisesYield(pars, R(p))
dfdσ = ForwardDiff.gradient(yield)
# Initializing variables
dσ = params[1:6]
dλ = params[end]
# Total strain
σ_new = σ_begin + dσ
σ_shifted = vec(σ_new - X)
# Creating wrapper for gradient
dΨdσ = dfdσ(σ_shifted)[1:6]
# Calculating plastic strain rate
dϵp = dλ * dΨdσ
# Updating material parameters
ϵp_new = ϵp + dϵp
dp = sqrt(2 / 3 * double_contraction(dϵp))
p_test = p + dp
# new hardening variables
R_ = R(p_test)
X_ = fX(ϵp_new)
# Evaluating equations
function_1 = dσ - C * (dϵ - dϵp)
function_2 = vonMisesYield(σ_new - X_, R_)
[vec(function_1); function_2]
end
"""
This is a novice implementation for double contraction, a=b:c
Parameters
----------
a: Array{Float64, 6}
Returns
-------
Float
"""
function double_contraction(a; b=a)
indexes = [1, 2, 3, 4, 5, 6, 4, 5, 6]
summation = 0
for i in indexes
summation += a[i]*b[i]
end
summation
end
"""
Function which calculates the stress. Also handles if any yielding happens
Parameters
----------
dϵ: Array{Float64, 6}
Strain rate vector in Voigt notation
Δt: Float
time increment
σ: Array{Float64, 6}
Last stress vector in Voigt notation
C: Array{Float64, (6, 6)}
Material tensor
p: Float
Accumulated plastic strain
X: Array{Float64, 6}
Kinematic hardening tensor
R: Function
Calculates isotropic hardening as a function of accumulated plastic strain
α:Float
Accumulated kinematic evolution variable
dα: Function
Calculates kinematic evolution variables rate as a function of accumulated plastic slip rate
fX: Function
Calculates the kinematic
Returns
-------
Tuple
returns following parameters: p, X, α, σ. See Parameters for definitions
"""
function calculate_stress(dϵ, σ, C, p, X, R, fX, ϵp)
# Test stress
σ_tria = σ + C * dϵ
# Calculating yield
yield = vonMisesYield(σ_tria - X, R(p))
if yield > 0
# Yielding happened
# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
initial_guess = [vec(σ_tria - σ); 0.1]
f(σ_) = G(σ_, dϵ, C, σ, p, X, R, fX, ϵp)
df = ForwardDiff.jacobian(f)
# Calculating root
result = nlsolve(not_in_place(f, df), initial_guess).zero
# Extracting values
σ += result[1:6]
dλ = result[end]
# Wrapper for gradient
yield_f(σ_) = vonMisesYield(σ_, R(p))
dfdσ_ = ForwardDiff.gradient(yield_f)
dΨdσ = dfdσ_(vec(σ)-X)
# Stress rate and plastic strain rate
dϵᵖ = dλ * dΨdσ
ϵp = ϵp + dϵᵖ
# Updating matrial parameters
dp = sqrt(2/3 * double_contraction(dϵᵖ))
p += dp
X = fX(ϵp)
else
σ = σ_tria
end
return (p, X, σ, ϵp)
end
Out[178]:
In [192]:
steps = 1000
strain_max = 0.003
num_cycles = 5
ϵ_tot = zeros(Float64, (steps, 6))
ϵ_tot2 = zeros(Float64, (steps, 6))
ϵ_tot3 = zeros(Float64, (steps, 6))
# Adding only strain in x-axis and counting for the poisson effect
ϵ_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
ϵ_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
ϵ_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
println("Done")
In [193]:
function isotropic_hardening(ϵp_cum, R0, Q)
return R0 + Q * ϵp_cum
end
function kinematic_hardening(ϵp, C)
return 2/3 * C * ϵp
end
Out[193]:
In [195]:
ϵ_last = zeros(Float64, (6)) # Last total strain
ϵᵖ = zeros(Float64, (6)) # Plastic strain
σ = zeros(Float64, (6, 1)) # Stress
ss = Float64[] # plotting stress
ee = Float64[] # plotting strain
pp = Float64[] # plotting strain
# Isotropic hardening
σy = 200.0 # yield limit
# Kinematic hardening
Ck = 2000.0 # saturation hardening C/D
X = zeros(Float64, 6) # kinematic hardening tensor
p = 0.0 # Accumulated plastic strain
# Isotropic hardening
R(ϵp_cum) = isotropic_hardening(ϵp_cum, σy, σy)
# kinematic hardening
fX(ϵp) = kinematic_hardening(ϵp, Ck)
for i=1:steps
# Actual calculation
dϵ = reshape(ϵ_tot[i, :, :], (6, 1)) - ϵ_last
p, X, σ, ϵᵖ = calculate_stress(dϵ, σ, C, p, X, R, fX, ϵᵖ)
ϵ_last += dϵ
push!(ss, σ[1])
push!(ee, ϵ_last[1])
push!(pp, p)
end
PyPlot.plot(ee, ss)
PyPlot.title("Stress-Strain curve")
PyPlot.xlabel("Strain")
PyPlot.ylabel("Stress")
PyPlot.grid()
In [ ]:
In [ ]: